Tutorial Part IV — Getting started with animal movement
Introduction
To start with, I would like to express my deep thanks to Cedric
Scherer for his help in designing and ameliorating the course scripts.
In this tutorial, we will show you how to get started with (real) animal movement data. To put it simple, animal movement data are just points in time and space, i.e. re-locations of an individual that can be provided as a data.frame with four columns. An identifier for the individual is only needed if you have various animals in one data set, i.e. if you have collared several individuals with the same collar or already appended several data sets:
| Identifier | Timestamp | x_coordinate | y_coordinate |
|---|---|---|---|
| Fox_1 | 2022-01-01 08:10:04 | 13.47027 | 52.497802 |
| Fox_1 | 2022-01-01 08:15:04 | 13.47015 | 52.497813 |
From Course2_SpatialR in our teaching
collection, you might remember that such data sets can be imported
as simple .txt or .csv files. Once you know in which coordinate
reference system (CRS) your coordinates are stemming from, you can
assign it to the data.frame, thereby you are
creating a geo-referenced data set: a
SpatialPointsDataframe-Object (if you work with the old
R-package {sp}) or an sf-Object (with package
{sf}). Looking at these x- and y-coordinates, you might
remember that they look like angular units, e.g. decimal degrees, and
the x-coordinate refers to longitude and the y-coordinate to latitude.
Hence, we are dealing with EPSG-code 4326. Once you have imported and
geo-referenced your data, you can plot and thoroughly check the
data.
Checking your data is a mandatory prerequisite
before any analysis (see Zuur
et al. 2010). With your movement data, you might have outliers in
space (e.g. when you test your collars in location A which might be
hundreds of kilometers from your actual trapping site), your data might
not have been regularly sampled because the animal was hiding and there
was no access to the satellite (weak GPS-signal) etc etc. This uneven
sampling can affect calculations of speed, for example, and therefore a
series of packages have been developed to deal with those issues in
movement data, like {move} or {amt}.
The most prominent movement analyses comprise home range estimation, calculations of habitat preference, behavioral analyses (e.g. Hertel et al. 2021) and detection of movement syndromes (personality differences) (Michelangeli et al. 2021). Some of the packages listed in the next chapter are designed for these analyses.
In the following, we will step by step start with loading and
exploring movement data of the female red fox (Vulpes vulpes)
‘Q von Stralau’, who had a small home range in Berlin, Germany. She was
collared in January 2018 and had a very stable daily routine as shown by
4 (20) min relocation intervals. Her location data are stored in the
file with the tag number tag5334_gps.txt. For any details
on the data please refer to Kimmig
(2021) and Scholz
(2020).
We will, however, only use the data of one month in this exercise, as
data sets quickly get too big. Please note: the courses Course1_IntroR
and Course2_SpatialR of our course script
collection
are obligatory for this tutorial.
Useful (web)sites and reading
- For analysis of telemetry data: packages adehabitatHR,
adehabitatLT, move, move2, recurse, momentuHMM, moveHMM, ctmm,
amt
Methods papers
- Joo, R, Boone, ME, Clay, TA, Patrick, SC, Clusella-Trullas, S, Basille, M. Navigating through the r packages for movement. J Anim Ecol. 2020; 89: 248– 267. https://doi.org/10.1111/1365-2656.13116
- Zuur, A.F., Ieno, E.N. and Elphick, C.S. (2010), A protocol for data exploration to avoid common statistical problems. Methods in Ecology and Evolution, 1: 3-14. https://doi.org/10.1111/j.2041-210X.2009.00001.x
Example papers
Hertel, AG, Royauté, R, Zedrosser, A, Mueller, T. Biologging reveals individual variation in behavioural predictability in the wild. J Anim Ecol. 2021; 90: 723– 737. https://doi.org/10.1111/1365-2656.13406
Kimmig, S (2021). The ecology of red foxes (Vulpes vulpes) in urban environments, PhD thesis, FU Berlin. https://refubium.fu-berlin.de/handle/fub188/32478
Michelangeli, M., Payne, E., Spiegel, O., Sinn, D. L., Leu, S. T., Gardner, M. G., & Sih, A. (2021). Personality, spatiotemporal ecological variation and resident/explorer movement syndromes in the sleepy lizard. Journal of Animal Ecology, 00, 1– 14. https://doi.org/10.1111/1365-2656.13616
Scholz, C (2020). The ecology of red foxes (Vulpes vulpes) in anthropogenic landscapes. PhD thesis, FU Berlin. https://refubium.fu-berlin.de/handle/fub188/29984
Getting started
To follow the tutorial, you can either clone or download the repository or you create your own R-project, copy the raw data and type the code chunks into an R-Script. Please refer to the section on using R-projects in Course2_RSpatial of our teaching collection.
If you start with your own R-project, I strongly recommend to use the
{d6}-package
Cedric Scherer provided. This package automatically sets up the ideal
folder structure: https://github.com/EcoDynIZW/d6
In any case, the course folder has the following structure:
#.
#└── d6_teaching_collection – (root folder)
# ├─── data
# ├─── data_move – (contains the file with the GPS locations)
# ├─── docs
# ├─── output – (contains the cleaned files and processed data)
# ├─── plots
# ├─── R – (contains the R-script)
# └ d6_teaching_collection.Rproj – (the RStudio-Project file location)
Necessary packages to install and load
We first have to install the packages and load them before we can use the functions we need.
## package names:
pkgs = c("here", "lubridate", "sf", "sp", "ggplot2", "tmap", "circular", "move",
"viridis", "devtools", "plotly","amt", "terra","effects","tidyterra",
"corrplot","ctmm","sjPlot","jtools")
# install.packages(pkgs) # only run this line once! for installing the packages!
# update.packages()Tip of the day: If you have already installed some of the packages above, you can first check which ones are already installed, and save the ones not installed in an object called ‘my_packages’ and only install the missing ones:
## character(0)
Now load the packages:
library(here) ## for easy directory management
library(lubridate) ## for date handling
library(sf) ## for handling simple feature objects
library(terra) ## for handling raster objects
library(ggplot2) ## for nice static plots
library(tmap) ## for nice interactive maps
library(circular) ## for circular stats + plots
library(move) ## for `move` objects
library(amt) ## for `move` objects
library(ctmm) ## for fitting continuous-time movement models
library(effects) ## for visualizing model results
#library(sjPlot) ## easy plotting of model results
library(jtools) ## easy plotting of model results
library(viridis) ## for perceptually uniform color palettes
library(corrplot) ## for correlation plotSet the working environment
Now that we are working inside an R-project, we can use the easy
functionality of the {here} package, which is automatically
pointing to your project’s root folder, e.g.:
Hence, there is no need to use the function setwd() any
more. Note: if it does not work, please close RStudio, go to your
Explorer and double-click on the .Rproj file. Then, under ‘files’
(usually lower right panel) double-click on the R folder and open the
script.
Load data
The movement data are stored in the data-raw subfolder. Let’s check which files are available.
## [1] ".../d6_teaching_collection/data/data_move/animal_relocations_32633.cpg"
## [2] ".../d6_teaching_collection/data/data_move/animal_relocations_32633.dbf"
## [3] ".../d6_teaching_collection/data/data_move/animal_relocations_32633.prj"
## [4] ".../d6_teaching_collection/data/data_move/animal_relocations_32633.qpj"
## [5] ".../d6_teaching_collection/data/data_move/animal_relocations_32633.shp"
## [6] ".../d6_teaching_collection/data/data_move/animal_relocations_32633.shx"
## [7] ".../d6_teaching_collection/data/data_move/geo-raw"
## [8] ".../projects_github/d6_teaching_collection/data/data_move/KettleHoles.txt"
## [9] ".../d6_teaching_collection/data/data_move/tag5334_gps.txt"
The output lists two results, there is another subfolder, and many
files. E.g., the element 7 of the vector lf, lf[7], is a
folder. If you already know that your movement data contains e.g. ‘gps’
in its name or is stored as ‘.txt’ or ‘.csv’ files, you can directly
search for those files with the pattern argument:
## check the difference, and note: full.names is set to FALSE
thefile <- list.files(path = here("data","data_move"), pattern = "gps", full.names = FALSE)
thefile## [1] "tag5334_gps.txt"
## ...and here to TRUE
thefullfile <- list.files(path = here("data","data_move"), pattern = "gps", full.names = TRUE)
thefullfile## [1] "C:/Users/admin/d6_teaching_collection/data/data_move/tag5334_gps.txt""
Now load the data file. Our first fox filename should be
tag5334_gps.txt, which we stored under object
‘thefullfilename’
Data check and cleaning
Checking for missing or incomplete information
Let’s have a look at the data. This is the typical way data are stored on e-obs collars:
## key.bin.checksum tag.serial.number start.timestamp longitude latitude
## 1 2135768252 5334 2018-10-30 04:00:00.000 13.47047 52.49499
## 2 2210577241 5334 2018-10-30 04:20:00.000 13.47157 52.49526
## 3 1824790132 5334 2018-10-30 04:40:00.000 13.47103 52.49525
## 4 4079427513 5334 2018-10-30 07:00:00.000 13.47109 52.49526
## 5 380027806 5334 2018-10-30 08:00:00.000 13.47105 52.49527
## height.above.ellipsoid type.of.fix status used.time.to.get.fix
## 1 76.3 3 A 22
## 2 81.6 3 A 18
## 3 66.5 3 A 53
## 4 73.5 3 A 24
## 5 71.1 3 A 47
## timestamp.of.fix battery.voltage fix.battery.voltage temperature
## 1 2018-10-30 04:00:23.000 3604 3389 10
## 2 2018-10-30 04:20:19.000 3616 3333 11
## 3 2018-10-30 04:40:54.000 3624 3374 13
## 4 2018-10-30 07:00:25.000 3659 3491 11
## 5 2018-10-30 08:00:48.000 3660 3470 13
## speed.over.ground heading.degree speed.accuracy.estimate
## 1 1.64 233.31 1.67
## 2 4.43 327.02 3.79
## 3 0.07 0.00 0.73
## 4 0.01 0.00 0.65
## 5 0.04 0.00 0.89
## horizontal.accuracy.estimate
## 1 7.94
## 2 59.65
## 3 15.36
## 4 7.17
## 5 8.45
For now, we will only work with few columns.
tag.serial.number refers to the individual
identifier (collar ID). There are two columns with timestamps:
start.timestamp is the preprogrammed
time-interval. Then there is the
timestamp.of.fix, which is the real time
the GPS-location was recorded. This is usually a bit later, i.e.
( =
start.timestamp +
used.time.to.get.fix + 1 second), as it
takes some time for the collar unit to connect to the satellite. The
spatial info is stored in the columns
longitude and
latitude.
Before you can transform the data.frame
dat_anim into a georeferenced spatial object, you need to
check whether there are missing locations in your data.frame,
otherwise you will get an error message on transformation.
This can happen if no GPS-signal could be recorded. Or -
importantly - some collars are only activated when the animal is moving
to save battery life. In that case, the missing GPS coordinates would
correspond with the last position (= be the same). Depending on which
analysis you want to do, e.g. define resting places, you might need to
fill the missing positions again.
In our case, there are a lot of missing values in the locations:
## if the latitude-entry is missing, the longitude value will also be missing
## so it is enough to only check the latitude
which(is.na(dat_anim$latitude))## [1] 48 65 89 90 91 129 183 185 222 224 226 229 260 266 355
## [16] 358 415 421 422 423 457 474 510 511 577 601 616 621 626 646
## [31] 673 674 685 696 721 774 784 796 822 833 842 855 863 907 921
## [46] 940 993 1006 1038 1070 1077 1078 1088 1126 1150 1159 1183 1188 1215 1238
## [61] 1268 1272 1304 1330 1480 1516 1550 1552 1559 1600 1637 1680 1717 1721 1736
## [76] 1738 1742 1772 1820 1824 1830 1867 1884 1973 1997 1999 2020 2023 2024 2025
## [91] 2026 2064 2075 2076 2109 2162 2234 2313 2329 2374 2379 2396 2419 2453 2455
## [106] 2518 2551 2609 2614 2638 2687 2728 2745 2770 2776 2783 2800 2804 2805 2854
## [121] 2926 2927 2931 2935 2943 2950 2961 2963 2968 2990 2991 3043 3044 3052 3090
## [136] 3102 3114 3172 3225 3226 3305 3328 3361 3372 3373 3386 3409 3498 3516 3517
## [151] 3518 3520 3591 3622 3623 3624
Delete rows with missing spatial info:
Make the crosscheck if there is missing info in a row in longitude. This should NOT be the case after we had deleted those rows:
## integer(0)
Checking for coarse spatial outliers
It could happen the collar was tested e.g. in Berlin, but the animal was finally caught and collared far away. Make a quick check whether there are strange locations:
There do not seem to be coarse outliers, data look compact.
Working with the date format
We will now add some additional columns, where separate days
(numbered from 1 to 365), the month (from 1 to 12) and the hour of the
day are stored (from 0 to 23). Sunset and sunrise can be calculated
based on dates and the location (latitude, longitude). This can be done
with the package {lubridate}. Check the vignette!
## have a look at the timestamp format, here the first row of the column start.timestamp
dat_anim_na$start.timestamp[1]## [1] "2018-10-30 04:00:00.000"
## define date-time format - the format is year-month-day_hour:min:sec:
dat_anim_na$start.timestamp <- ymd_hms(dat_anim_na$start.timestamp, tz="Europe/Berlin")
dat_anim_na$start.timestamp[1]## [1] "2018-10-30 04:00:00 CET"
Now append the information:
dat_anim_na$yearday <- yday(dat_anim_na$start.timestamp)
dat_anim_na$month <- month(dat_anim_na$start.timestamp)
dat_anim_na$hour <- hour(dat_anim_na$start.timestamp)
dat_anim_na$kweek <- week(dat_anim_na$start.timestamp)
dat_anim_na$date <- date(dat_anim_na$start.timestamp)
## crosscheck with
# head(dat_anim_na)In addition, we want to calculate the hours of sunset and sunrise as well as daylength. For this, we need to install a package that is still under development, i.e. which is not on CRAN. We therefore must download and install it locally:
# devtools::install_github("bgctw/solartime") ## run only once to install package from GitHub
library(solartime)For computing sunset and sunrise, the latitude and longitude must be provided as well:
dat_anim_na$sunrise <- computeSunriseHour(timestamp = dat_anim_na$start.timestamp,
latDeg = dat_anim_na$latitude,
longDeg = dat_anim_na$longitude)
dat_anim_na$sunset <- computeSunsetHour(dat_anim_na$start.timestamp,
dat_anim_na$latitude,
dat_anim_na$longitude)
dat_anim_na$daylength <- computeDayLength(dat_anim_na$start.timestamp,dat_anim_na$latitude)
dat_anim_na$daytime <- computeIsDayByLocation(dat_anim_na$start.timestamp,
dat_anim_na$latitude,
dat_anim_na$longitude)
## check new variables
# head(dat_anim_na)
hist(dat_anim_na$daylength)unique(dat_anim_na$daytime) ## gives the levels, i.e. 'TRUE' means daylight, 'FALSE' means nighttime and darkness## [1] FALSE TRUE
Check for temporal outliers
Check if there are strange dates, or dates before you collared the animal (e.g., usually the collar is activated and tested before the animal is collared):
##
## 2018-10-30 2018-10-31 2018-11-01 2018-11-02 2018-11-03 2018-11-04 2018-11-05
## 24 40 35 39 42 34 34
## 2018-11-06 2018-11-07 2018-11-08 2018-11-09 2018-11-10 2018-11-11 2018-11-12
## 41 42 42 41 38 38 39
## 2018-11-13 2018-11-14 2018-11-15 2018-11-16 2018-11-17 2018-11-18 2018-11-19
## 41 38 34 38 44 34 35
## 2018-11-20 2018-11-21 2018-11-22 2018-11-23 2018-11-24 2018-11-25 2018-11-26
## 36 34 40 46 40 37 41
## 2018-11-27 2018-11-28 2018-11-29 2018-11-30 2018-12-01 2018-12-02 2018-12-03
## 43 39 34 42 45 37 39
## 2018-12-04 2018-12-05 2018-12-06 2018-12-07 2018-12-08 2018-12-09 2018-12-10
## 38 39 38 31 38 36 36
## 2018-12-11 2018-12-12 2018-12-13 2018-12-14 2018-12-15 2018-12-16 2018-12-17
## 39 34 38 37 37 43 43
## 2018-12-18 2018-12-19 2018-12-20 2018-12-21 2018-12-22 2018-12-23 2018-12-24
## 37 31 41 36 38 42 44
## 2018-12-25 2018-12-26 2018-12-27 2018-12-28 2018-12-29 2018-12-30 2018-12-31
## 32 30 38 38 36 37 39
## 2019-01-01 2019-01-02 2019-01-03 2019-01-04 2019-01-05 2019-01-06 2019-01-07
## 33 43 37 36 39 39 45
## 2019-01-08 2019-01-09 2019-01-10 2019-01-11 2019-01-12 2019-01-13 2019-01-14
## 36 38 41 33 26 35 40
## 2019-01-15 2019-01-16 2019-01-17 2019-01-18 2019-01-19 2019-01-20 2019-01-21
## 28 35 35 40 39 34 31
## 2019-01-22 2019-01-23 2019-01-24 2019-01-25 2019-01-26 2019-01-27 2019-01-28
## 36 36 34 36 43 34 33
## 2019-01-29 2019-01-30 2019-01-31 2025-12-26
## 35 34 27 1
It might be easier to spot visually. We use the
{ggplot2} package here as it recognizes and respects the
date format:
There is a strange date → 2025-12-26!
Delete this data row:
## key.bin.checksum tag.serial.number start.timestamp longitude latitude
## 2622 2948730959 5334 2025-12-26 04:00:00 13.47851 52.49133
## height.above.ellipsoid type.of.fix status used.time.to.get.fix
## 2622 78.8 3 A 23
## timestamp.of.fix battery.voltage fix.battery.voltage temperature
## 2622 2025-12-26 04:00:24.000 3582 3287 3
## speed.over.ground heading.degree speed.accuracy.estimate
## 2622 0.03 0 0.64
## horizontal.accuracy.estimate yearday month hour kweek date sunrise
## 2622 4.35 360 12 4 52 2025-12-26 8.390617
## sunset daylength daytime
## 2622 15.80777 7.417151 FALSE
dat_anim_na <- dat_anim_na[-delme,] ## delete the strange date and
table(dat_anim_na$date) ## check again##
## 2018-10-30 2018-10-31 2018-11-01 2018-11-02 2018-11-03 2018-11-04 2018-11-05
## 24 40 35 39 42 34 34
## 2018-11-06 2018-11-07 2018-11-08 2018-11-09 2018-11-10 2018-11-11 2018-11-12
## 41 42 42 41 38 38 39
## 2018-11-13 2018-11-14 2018-11-15 2018-11-16 2018-11-17 2018-11-18 2018-11-19
## 41 38 34 38 44 34 35
## 2018-11-20 2018-11-21 2018-11-22 2018-11-23 2018-11-24 2018-11-25 2018-11-26
## 36 34 40 46 40 37 41
## 2018-11-27 2018-11-28 2018-11-29 2018-11-30 2018-12-01 2018-12-02 2018-12-03
## 43 39 34 42 45 37 39
## 2018-12-04 2018-12-05 2018-12-06 2018-12-07 2018-12-08 2018-12-09 2018-12-10
## 38 39 38 31 38 36 36
## 2018-12-11 2018-12-12 2018-12-13 2018-12-14 2018-12-15 2018-12-16 2018-12-17
## 39 34 38 37 37 43 43
## 2018-12-18 2018-12-19 2018-12-20 2018-12-21 2018-12-22 2018-12-23 2018-12-24
## 37 31 41 36 38 42 44
## 2018-12-25 2018-12-26 2018-12-27 2018-12-28 2018-12-29 2018-12-30 2018-12-31
## 32 30 38 38 36 37 39
## 2019-01-01 2019-01-02 2019-01-03 2019-01-04 2019-01-05 2019-01-06 2019-01-07
## 33 43 37 36 39 39 45
## 2019-01-08 2019-01-09 2019-01-10 2019-01-11 2019-01-12 2019-01-13 2019-01-14
## 36 38 41 33 26 35 40
## 2019-01-15 2019-01-16 2019-01-17 2019-01-18 2019-01-19 2019-01-20 2019-01-21
## 28 35 35 40 39 34 31
## 2019-01-22 2019-01-23 2019-01-24 2019-01-25 2019-01-26 2019-01-27 2019-01-28
## 36 36 34 36 43 34 33
## 2019-01-29 2019-01-30 2019-01-31
## 35 34 27
Save the cleaned file
Finally, we save the processed data file that we will use for
exploration and analysis into the subfolder
.../output/data-proc. There are two options we can use:
- R data file - can only be opened/ read with R by using function readRDS()
- interchange file format .csv
Check your output folder for these files. The efficient ‘.Rds’ file
storage is about 3 times smaller than the ‘.csv’ file.
Exercise 4.1
This is a recap from Course2. Please plot the animal
relocations in two different colours based on the column
daytime, using one of the options to create maps with
{ggplot2}, {ggmap}, {leaflet} ,
{tmap} or {ggspatial}. Use the processed file
of fox Q (tag5334_gps_proc) and do it in a separate script.
Save your script as Course4_Exercise1_*yourname*.R.
Hint: Remember to load the relevant libraries. Note that you
might want to plot the locations in the correct spatial dimensions by
projecting it using the functions st_as_sf() and
st_transform().
Data exploration
Load the cleaned data file
I recommend that you store the raw file safely and continue with the cleaned and processed file after major data manipulations have been conducted to minimize errors. I’d even suggest to make these steps in different R-scripts.
Have a look whether there are gaps/ missing days in the data, or
whether we have approximately a regular number of fixes each day. We can
use {ggplot2} to plot true dates:
Now let’s plot the Julian day, yearday:
The number of fixes (locations taken) per day seems to be quite regular. Mind the difference when plotting per date (sorted) and yearday (1-365), irrespective of the year.
Plot activity
Now we can have a look at the distribution of the fixes during the day, which is a hint on the active phase of the animal:
More fixes during the night - this is the active phase. As these are circular data, we can plot the number of fixes as sign for activity (knowing that our collars did not record when foxes were inactive! Otherwise, we cannot distinguish missing data from inactive times! -> bias).
## simple circular plot (rose diagram)
timetoplot <- circular(anim_proc$hour %% 24, ## convert to 24 hrs = bins
units = "hours", template = "clock24")
## Note: use namespace `circular::` here as there are multiple functions called `rose.diag()`
circular::rose.diag(timetoplot, bin = 24, col = "blue",
main = "Events by Hour (sqrt scale)", prop = 3)Check whether the activity during the day was before sunrise (our ‘daytime’ column, which is either TRUE (=yes, daylight) or FALSE (dark hours):
# with ggplot
# code adapted from https://gist.github.com/mattbaggott/4361381
ggplot(anim_proc, aes(x = hour,fill = daytime)) +
geom_bar(width = 1, color = "grey20") +
coord_polar(start = -0.15) + ## to center setchange start
theme_minimal() +
scale_x_continuous(breaks = 0:23) +
scale_fill_brewer(palette = "Set2", name = NULL, labels = c("Night", "Day")) +
labs(x = NULL, y = "Count", title = "Events by Time of the Day")Apart from the outliers between 11 and 13 o’clock, the active phase was during the dark hours.
Actogramm
Convert the data into a spatial object
Now we transform the data into sf and sp
objects and assign the coordinate reference system:
# transform into spatial simple feature sf object
mydf_sf <- st_as_sf(x = data.frame(anim_proc),
coords = c("longitude", "latitude"),
crs = 4326,
sf_column_name = "geometry" )
# transform into SpatialPointsDataFrame - for crosschecking
mydf_sp <- as(mydf_sf, "Spatial") And then we project the reference system from angular units to a planar coordinate reference system in meters:
# transform CRS to projected one in meter distance units
mydf_sf_trans <- st_transform(mydf_sf, 3035 ) # EPSG-code
mydf_sp_trans <- spTransform(mydf_sp, CRS("+init=epsg:3035")) ## TODO warning ##skip spRecently, there are issues with the missing datum in the
CRS-specifications. We ignore this for now. More on the issue of moving
from proj4 to proj6 in the future:
* https://inbo.github.io/tutorials/tutorials/spatial_crs_coding/
* https://cran.r-project.org/web/packages/sp/vignettes/CRS_warnings.html
Make a quick plot of the data
Recap for styling the plot:
https://www.r-bloggers.com/2021/12/introduction-to-geospatial-visualization-with-the-tmap-package/
tmap_mode(mode = "view")
tm_shape(shp = mydf_sf_trans) +
tm_dots(size = 0.01,
col = "daytime",
alpha = 0.5)Have a deep look at the data: did the fox really swim? (points in
Lake Rummelsburg). Or are these outliers? No: the lake was frozen -
check the date! Note: you really need to know your animals and the area
before analysing data!
Basic movement metrics
Transform your data into an {amt} object
Check here for the possibilities of what to do with the
{amt}-package:
* https://cran.r-project.org/web/packages/amt/vignettes/p1_getting_started.html
Because we are dealing with spatial measures (e.g. area, distance, perimeter,…), I recommend to work with projected coordinates. To create a move-object, the coordinates of the projected CRS need to be appended to a data.frame:
## assign columns with coordinates_epsg to add later into the data frame for amt-object
mydf_sf_trans$X_3035 <- st_coordinates(mydf_sf_trans)[,1]
mydf_sf_trans$Y_3035 <- st_coordinates(mydf_sf_trans)[,2]
## new object as data frame
anim_for_move <- st_drop_geometry(mydf_sf_trans)
head(anim_for_move)## key.bin.checksum tag.serial.number start.timestamp height.above.ellipsoid
## 1 2135768252 5334 2018-10-30 04:00:00 76.3
## 2 2210577241 5334 2018-10-30 04:20:00 81.6
## 3 1824790132 5334 2018-10-30 04:40:00 66.5
## 4 4079427513 5334 2018-10-30 07:00:00 73.5
## 5 380027806 5334 2018-10-30 08:00:00 71.1
## 6 3352806116 5334 2018-10-30 12:00:00 71.5
## type.of.fix status used.time.to.get.fix timestamp.of.fix
## 1 3 A 22 2018-10-30 04:00:23.000
## 2 3 A 18 2018-10-30 04:20:19.000
## 3 3 A 53 2018-10-30 04:40:54.000
## 4 3 A 24 2018-10-30 07:00:25.000
## 5 3 A 47 2018-10-30 08:00:48.000
## 6 3 A 23 2018-10-30 12:00:24.000
## battery.voltage fix.battery.voltage temperature speed.over.ground
## 1 3604 3389 10 1.64
## 2 3616 3333 11 4.43
## 3 3624 3374 13 0.07
## 4 3659 3491 11 0.01
## 5 3660 3470 13 0.04
## 6 3664 3508 14 0.02
## heading.degree speed.accuracy.estimate horizontal.accuracy.estimate yearday
## 1 233.31 1.67 7.94 303
## 2 327.02 3.79 59.65 303
## 3 0.00 0.73 15.36 303
## 4 0.00 0.65 7.17 303
## 5 0.00 0.89 8.45 303
## 6 0.00 0.44 3.84 303
## month hour kweek date sunrise sunset daylength daytime X_3035
## 1 10 4 44 2018-10-30 7.033396 16.62727 9.593878 FALSE 4556608
## 2 10 4 44 2018-10-30 7.033335 16.62719 9.593855 FALSE 4556681
## 3 10 4 44 2018-10-30 7.033371 16.62723 9.593855 FALSE 4556645
## 4 10 7 44 2018-10-30 7.033367 16.62722 9.593855 FALSE 4556649
## 5 10 8 44 2018-10-30 7.033370 16.62722 9.593853 TRUE 4556646
## 6 10 12 44 2018-10-30 7.033371 16.62722 9.593852 TRUE 4556646
## Y_3035
## 1 3270698
## 2 3270731
## 3 3270729
## 4 3270729
## 5 3270731
## 6 3270732
## [1] TRUE
my_fox <- make_track(tbl = anim_for_move,
.x = X_3035,
.y = Y_3035,
.t = start.timestamp,
id = tag.serial.number, ## has to be a column of data set
crs = 3035) ## has to be set as EPSG code
# all_cols = TRUE) ## if all columns need to be kept
head(my_fox)## # A tibble: 6 × 4
## x_ y_ t_ id
## * <dbl> <dbl> <dttm> <int>
## 1 4556608. 3270698. 2018-10-30 04:00:00 5334
## 2 4556681. 3270731. 2018-10-30 04:20:00 5334
## 3 4556645. 3270729. 2018-10-30 04:40:00 5334
## 4 4556649. 3270729. 2018-10-30 07:00:00 5334
## 5 4556646. 3270731. 2018-10-30 08:00:00 5334
## 6 4556646. 3270732. 2018-10-30 12:00:00 5334
## [1] FALSE
A hint (for the exercises later, so please read!): check the argument ‘id’. Here, I put one tag-id number, however, if multiple individuals are used, you can provide a column with multiple IDs.
Movement metrics
## Time lags between the fixes:
summarize_sampling_rate(my_fox) ## varies a lot, but median here 20 min, max 4 hrs## # A tibble: 1 × 9
## min q1 median mean q3 max sd n unit
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <int> <chr>
## 1 17.2 20 20 38.5 20 240 46.8 3502 min
## The distance ranges form 0.5 m to 1 km,
## but we have here the 4 hrs time lags included
range(step_lengths(my_fox), na.rm =TRUE)## [1] 0.5468472 1028.7766731
## [1] "track_xyt" "track_xy" "tbl_df" "tbl" "data.frame"
#max(amt::speed(my_fox),na.rm=TRUE) ## 0.8, measured in m/s
## Let's look only at 60 min intervals:
my_fox_60 <- track_resample(my_fox, rate = minutes(60), tolerance = minutes(5))
range(step_lengths(my_fox_60), na.rm =TRUE)## [1] 0.7395256 1184.1579983
## can you explain why the step length gets larger?
## what happens with speed?
#max(speed(my_fox_60),na.rm=TRUE) Turning angles:
The calculation of turning angles only makes
sense with very high resolution data. With 20 min intervals between
fixes, turning angles do not make much sense any more, as the animals
can have turned and moved a lot in between. I nevertheless show the
code, because turning angles are often used to describe behavior.
Signer et al. 2018; DOI: 10.1002/ece3.4823
We will also choose to
keep only those bursts (subsets of the track with constant sampling
rate, within the specified tolerance) with at least three relocations,
the minimum required to calculate a turn angle (amt::filter_min _ n _
burst). The following code implements those choices and translates from
a point representation to a step (step length, turn angle)
representation of the data.
## first transform track into steps
tmp1 <- amt::filter_min_n_burst(my_fox_60, min_n=3)
## then derive turning angles - column 'ta_' will be appended
fox_steps_60 <- steps_by_burst(tmp1) ## turn. ang. from -pi to pi = -180° to 180°
range(as_degree(fox_steps_60$ta_), na.rm = TRUE) ## convert to degree## [1] -179.9449 179.9833
The homerange concept - MCP, kernelUD (aKDE)
See lecture!
Minimum Convex Polygon
Get the home range size.
https://cran.r-project.org/web/packages/amt/vignettes/p2_hr.html
mcp_fox_95 <- amt::hr_mcp(x = my_fox, ## the track object
levels = c(0.95)) # define the level of the mcp.
plot(mcp_fox_95) ## quick and uggly## if you are not sure about how to get the info of the object, use str()
mcp_fox_95$mcp$area/ 1e6 ## area converted to km2## 0.586523 [m^2]
## # A tibble: 1 × 3
## level what area
## <dbl> <chr> <dbl>
## 1 0.95 estimate 586523.
## calculated MCP at different levels
mcp_fox <- hr_mcp(my_fox, levels = seq(0.2, 1, 0.1))
hr_area(mcp_fox)## # A tibble: 9 × 3
## level what area
## <dbl> <chr> <dbl>
## 1 1 estimate 684232.
## 2 0.9 estimate 509770.
## 3 0.8 estimate 305613.
## 4 0.7 estimate 232182.
## 5 0.6 estimate 162707.
## 6 0.5 estimate 131338.
## 7 0.4 estimate 123252.
## 8 0.3 estimate 110913.
## 9 0.2 estimate 79550.
Convert to sf object to store as shapefile, for example, and make a
plot
mcp_fox_95_sf <- sf::st_as_sf(mcp_fox_95$mcp)
st_write(mcp_fox_95_sf,
dsn = here("output", "mcp_fox_95_sf.shp"),
delete_layer = TRUE)## Deleting layer `mcp_fox_95_sf' using driver `ESRI Shapefile'
## Writing layer `mcp_fox_95_sf' to data source
## `C:/Users/kramer/_github_lokal_NB0144/d6_teaching_collection/output/mcp_fox_95_sf.shp' using driver `ESRI Shapefile'
## Writing 1 features with 3 fields and geometry type Polygon.
ggplot() +
geom_sf(data = mcp_fox_95_sf) +
geom_sf(data = mydf_sf_trans, size=0.1, col = 'red') + ## df converted to sf above
coord_sf()Kernel utility density
Calculate the density kernel:
## # A tibble: 1 × 3
## level what area
## <dbl> <chr> <dbl>
## 1 0.95 estimate 532654.
## # A tibble: 1 × 3
## level what area
## <dbl> <chr> <dbl>
## 1 0.95 estimate 586523.
## spatial density distribution
# plot(kde_fox_95$ud,col =viridis(100, direction = -1)) ## density distribution
# contour(kde_fox_95$ud)Also here, convert to sf object to store as shapefile, for example,
and make a plot
kde_fox_95_sf <- hr_isopleths(kde_fox_95) ## nicely, this is already an sf object
st_write(kde_fox_95_sf,
dsn = here("output", "kde_fox_95_sf.shp"),
delete_layer = TRUE)## Deleting layer `kde_fox_95_sf' using driver `ESRI Shapefile'
## Writing layer `kde_fox_95_sf' to data source
## `C:/Users/kramer/_github_lokal_NB0144/d6_teaching_collection/output/kde_fox_95_sf.shp' using driver `ESRI Shapefile'
## Writing 1 features with 3 fields and geometry type Multi Polygon.
ggplot() +
geom_sf(data = mcp_fox_95_sf, alpha=0.5) +
geom_sf(data = kde_fox_95_sf, alpha=0.1) +
geom_sf(data = mydf_sf_trans, size=0.1, col = 'red') + ## df converted to sf above
coord_sf()Resource selection
Nothing makes sense but in the light of environmental
information…
https://cran.r-project.org/web/packages/amt/vignettes/p3_rsf.html
Let’s first thin and filter the dataset to one location per day
## Let's look only at 60 min intervals:
my_fox_24 <- track_resample(my_fox, rate = hours(24), tolerance = minutes(25))
dim(my_fox_24)## [1] 95 5
Environmental covariates
First, load some environmental data. Check what is available in your data_berlin-course folder, check for the CRS and assign if missing:
maps_wd <- here("data","data_berlin","geo_raster_current_gtif")
lf <- list.files(maps_wd, full.names = TRUE); lf## [1] "C:/Users/kramer/_github_lokal_NB0144/d6_teaching_collection/data/data_berlin/geo_raster_current_gtif/imperviousness_2012_100m_3035.tif"
## [2] "C:/Users/kramer/_github_lokal_NB0144/d6_teaching_collection/data/data_berlin/geo_raster_current_gtif/noise_daynight_2017_100m_3035.tif"
## [3] "C:/Users/kramer/_github_lokal_NB0144/d6_teaching_collection/data/data_berlin/geo_raster_current_gtif/summer_temp_100m_3035.tif"
## [4] "C:/Users/kramer/_github_lokal_NB0144/d6_teaching_collection/data/data_berlin/geo_raster_current_gtif/tree_cover_density_2012_100m_3035.tif"
## [5] "C:/Users/kramer/_github_lokal_NB0144/d6_teaching_collection/data/data_berlin/geo_raster_current_gtif/water_bodies_2010_100m_3035.tif"
b_imperv <- terra::rast(lf[1])
b_sutemp <- terra::rast(lf[3])
b_forest <- terra::rast(lf[4])
b_noise <- terra::rast(lf[2])
water_raster <- rast(here::here("data","data_berlin","geo_raster_current_gtif","water_bodies_2010_100m_3035.tif"))
## chech CRS:
b_sutemp ## coord. ref. is assigned because it is a geoTiff, not ascii-file.## class : SpatRaster
## dimensions : 570, 657, 1 (nrow, ncol, nlyr)
## resolution : 100, 100 (x, y)
## extent : 4521040, 4586740, 3243800, 3300800 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs
## source : summer_temp_100m_3035.tif
## name : summer_temp_100m_3035
## min value : 2052
## max value : 3133
## combine to raster stack
b_env_var_stack <- terra::rast(list(b_imperv,b_sutemp, b_forest,b_noise))
names(b_env_var_stack) <- c("imp","stmp","tree","noise")
terra::plot(b_env_var_stack) ## somehow does not work any more - TODOCompute a correlation matrix before you enter all variables into the model to check for multicollinearity
M <- terra::layerCor(b_env_var_stack, fun = "pearson")
corrplot(M$correlation, type="upper",tl.col="black", tl.srt=45)Now create random points in MCP home range of the fox. Default is 10 times more random points than points. It automatically sets a column ‘case_’ = FALSE, because it is a random point, not a true animal location.
Now we extract the covariates (several ways to do so, see Course 2 Spatial R and function ‘st_intersection()’ or ‘extract’, but we use the function ‘extract_covariates()’ of the amt-package here)
env_rp_fox_24 <- amt::extract_covariates(x=rp_fox_24, covariates=b_env_var_stack)
head(env_rp_fox_24)## # A tibble: 6 × 7
## case_ x_ y_ imp stmp tree noise
## * <lgl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 FALSE 4557031. 3270582. 9304 2674 4800 5236
## 2 FALSE 4557041. 3270513. 7078 2664 6400 5194
## 3 FALSE 4557321. 3270166. 3400 2683 7725 5252
## 4 FALSE 4556747. 3270704. 8620 2661 2420 4717
## 5 FALSE 4556626. 3270690. 9356 2693 3627 5321
## 6 FALSE 4556651. 3270681. 7783 2685 5892 5069
We do the same with the fox locations. Note that we have a different structure of the data.frame with more columns. And we have to add the column ‘case_’
## # A tibble: 6 × 9
## x_ y_ t_ id burst_ imp stmp tree noise
## * <dbl> <dbl> <dttm> <int> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 4556608. 3270698. 2018-10-30 04:00:00 5334 1 9356 2693 3627 5321
## 2 4556644. 3270734. 2018-10-31 03:40:00 5334 1 8753 2666 5271 5309
## 3 4556808. 3270550. 2018-11-01 03:40:00 5334 1 7175 2694 5650 5030
## 4 4557182. 3270347. 2018-11-02 03:20:00 5334 1 2150 2662 7216 5257
## 5 4557442. 3270296. 2018-11-03 03:00:00 5334 1 10000 2626 0 5035
## 6 4556524. 3270576. 2018-11-04 02:40:00 5334 1 8850 2726 0 5742
## [1] "x_" "y_" "t_" "id" "burst_" "imp" "stmp" "tree"
## [9] "noise" "case_"
## # A tibble: 6 × 7
## case_ x_ y_ imp stmp tree noise
## * <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 TRUE 4556608. 3270698. 9356 2693 3627 5321
## 2 TRUE 4556644. 3270734. 8753 2666 5271 5309
## 3 TRUE 4556808. 3270550. 7175 2694 5650 5030
## 4 TRUE 4557182. 3270347. 2150 2662 7216 5257
## 5 TRUE 4557442. 3270296. 10000 2626 0 5035
## 6 TRUE 4556524. 3270576. 8850 2726 0 5742
Combine datasets for analysis and convert binary TRUE/FALSE (column ‘case_’) into numeric.
## # A tibble: 6 × 7
## case_ x_ y_ imp stmp tree noise
## * <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 FALSE 4557031. 3270582. 9304 2674 4800 5236
## 2 FALSE 4557041. 3270513. 7078 2664 6400 5194
## 3 FALSE 4557321. 3270166. 3400 2683 7725 5252
## 4 FALSE 4556747. 3270704. 8620 2661 2420 4717
## 5 FALSE 4556626. 3270690. 9356 2693 3627 5321
## 6 FALSE 4556651. 3270681. 7783 2685 5892 5069
RSF model fit
Scaling variables, so that they all have the same range of values, can be helpful when plotting - only few plot functions below have inbuilt functions to standardize while plotting.
rsf_fox <- glm(data = df_rsf_fox_24,
formula = case_num ~ scale(noise) + scale(imp) + scale(stmp) + scale(tree),
#formula = case_num ~ noise + imp + stmp + tree,
family = binomial(link = "logit"))
summary(rsf_fox)##
## Call:
## glm(formula = case_num ~ scale(noise) + scale(imp) + scale(stmp) +
## scale(tree), family = binomial(link = "logit"), data = df_rsf_fox_24)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.73927 0.13540 5.460 4.76e-08 ***
## scale(noise) 0.40285 0.14868 2.710 0.006737 **
## scale(imp) 0.83588 0.17084 4.893 9.94e-07 ***
## scale(stmp) 0.04421 0.14000 0.316 0.752159
## scale(tree) 0.64029 0.17761 3.605 0.000312 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 362.81 on 284 degrees of freedom
## Residual deviance: 325.70 on 280 degrees of freedom
## AIC: 335.7
##
## Number of Fisher Scoring iterations: 4
Effects plots
Plot the estimates as effects plots and forest plots.
Forest plot
https://www.jtools.jacob-long.com/articles/summ.html
The r-package `{flextable}´ provides an easy way to export the model
summary from R into an excel or word file format.
To sum up: The
results - that a fox prefers houses and noise - do not really make
sense. What could be the reason be?
You could repeat the analysis
with day and nighttime locations of the fox….
Predict the model
## class : SpatRaster
## dimensions : 570, 657, 1 (nrow, ncol, nlyr)
## resolution : 100, 100 (x, y)
## extent : 4521040, 4586740, 3243800, 3300800 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs
## source(s) : memory
## name : lyr1
## min value : -15.96097
## max value : 10.51901
## TODO - still in logit scale? oder liegts am scale?
# plot map
terra::plot(fox_pred, col = viridis(100) )
terra::plot(water_raster, col = "darkslategray1", legend=FALSE, add = TRUE)-> go to Exercise 4.2;
check argument ‘nest: -id’ in
amt-package to work with several animals
Coming soon: Recurse analysis, HMM, …
END
Adding the session info can be very helpful when going back to old scripts or using scripts of others:
Session Info
## R version 4.4.2 (2024-10-31 ucrt)
## Platform: x86_64-w64-mingw32/x64
## Running under: Windows 10 x64 (build 19044)
##
## Matrix products: default
##
##
## locale:
## [1] LC_COLLATE=English_United Kingdom.utf8
## [2] LC_CTYPE=English_United Kingdom.utf8
## [3] LC_MONETARY=English_United Kingdom.utf8
## [4] LC_NUMERIC=C
## [5] LC_TIME=C
##
## time zone: Europe/Berlin
## tzcode source: internal
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] solartime_0.0.4 corrplot_0.95 viridis_0.6.5 viridisLite_0.4.2
## [5] jtools_2.3.0 effects_4.2-2 carData_3.0-5 ctmm_1.2.0
## [9] amt_0.2.2.0 move_4.2.6 raster_3.6-30 sp_2.1-4
## [13] geosphere_1.5-20 circular_0.5-1 tmap_3.3-4 ggplot2_3.5.1
## [17] terra_1.8-5 sf_1.0-19 lubridate_1.9.4 here_1.0.1
##
## loaded via a namespace (and not attached):
## [1] Rdpack_2.6.2 DBI_1.2.3 gridExtra_2.3
## [4] tmaptools_3.1-1 s2_1.1.7 rlang_1.1.4
## [7] magrittr_2.0.3 furrr_0.3.1 e1071_1.7-16
## [10] compiler_4.4.2 png_0.1-8 systemfonts_1.1.0
## [13] vctrs_0.6.5 wk_0.9.4 pkgconfig_2.0.3
## [16] fastmap_1.2.0 backports_1.5.0 labeling_0.4.3
## [19] lwgeom_0.2-14 rmdformats_1.0.4 leafem_0.2.3
## [22] utf8_1.2.4 pander_0.6.6 rmarkdown_2.29
## [25] nloptr_2.1.1 ragg_1.3.3 purrr_1.0.2
## [28] xfun_0.49 cachem_1.1.0 jsonlite_1.9.1
## [31] broom_1.0.7 parallel_4.4.2 R6_2.6.1
## [34] bslib_0.8.0 RColorBrewer_1.1-3 parallelly_1.42.0
## [37] boot_1.3-31 estimability_1.5.1 jquerylib_0.1.4
## [40] stars_0.6-7 Rcpp_1.0.14 bookdown_0.41
## [43] knitr_1.49 base64enc_0.1-3 leaflet.providers_2.0.0
## [46] Matrix_1.7-1 splines_4.4.2 nnet_7.3-19
## [49] timechange_0.3.0 tidyselect_1.2.1 rstudioapi_0.17.1
## [52] dichromat_2.0-0.1 abind_1.4-8 yaml_2.3.10
## [55] codetools_0.2-20 listenv_0.9.1 lattice_0.22-6
## [58] tibble_3.2.1 leafsync_0.1.0 withr_3.0.2
## [61] evaluate_1.0.1 future_1.34.0 survival_3.7-0
## [64] units_0.8-6 proxy_0.4-27 survey_4.4-2
## [67] xml2_1.3.6 pillar_1.10.0 KernSmooth_2.23-24
## [70] checkmate_2.3.2 insight_1.1.0 generics_0.1.3
## [73] rprojroot_2.0.4 munsell_0.5.1 scales_1.3.0
## [76] minqa_1.2.8 globals_0.16.3 class_7.3-22
## [79] glue_1.8.0 tools_4.4.2 data.table_1.16.4
## [82] lme4_1.1-35.5 forcats_1.0.0 mvtnorm_1.3-2
## [85] XML_3.99-0.17 grid_4.4.2 tidyr_1.3.1
## [88] mitools_2.4 rbibutils_2.3 crosstalk_1.2.1
## [91] colorspace_2.1-1 nlme_3.1-166 cli_3.6.3
## [94] textshaping_0.4.1 dplyr_1.1.4 gtable_0.3.6
## [97] broom.mixed_0.2.9.6 sass_0.4.9 digest_0.6.37
## [100] classInt_0.4-11 farver_2.1.2 htmlwidgets_1.6.4
## [103] memoise_2.0.1 htmltools_0.5.8.1 lifecycle_1.0.4
## [106] leaflet_2.2.2 httr_1.4.7 MASS_7.3-61